Intro
In this post I will map boundaries related to Memphis, TN.
Data Source
Boundaries are made available by the US Census Bureau on their TIGER/Line Shapefiles page. These files do not contain demographic data, but their GEOIDs can be linked to other Census data. In R, the tigris package can be used to access these boundaries.
I’ll also use the sf package to ensure consistent coordinates and to filter geographic areas. I use ggplot in this post to make maps.
If you’re following along in R, load these packages to get started.
Determine which boundaries are best for your project. A full list of TIGER/Line boundaries available through tigris is available here.
CRS
In essence, a coordinate reference system (crs) tells a map how to look. If your map doesn’t look right, like it’s skewed or warped or whatever, you probably need to reangle the map by setting the crs. If you plan to combine maps, defining the crs ensures projections are consistent.
If you do not define a crs, tigris and sf will default to 4269 (NAD 1983). Set the crs of a dataset using st_transform(crs = ####). You can check a table’s crs using st_crs(). For more info see the section on crs in Analyzing US Census Data.
The package crsuggest can help find the correct crs for your map. I also found the website epsg.io useful; for example, you can quickly see the boundaries for the crs code 6510 at https://epsg.io/6510. This crs was the top suggestion for my maps, but the projection excluded Tennessee and I found the difference minimal from the default. For simplicity, I stuck with the default.
Memphis Metro Area
To understand a place, we often need to look further than the city limits and also study the surrounding area. One way to do this is to analyze Census metro and micro areas.
A metropolitan area is centered around an urban core and contains 50,000 or more people. A micropolitan area contains at least 10,000, but less than 50,000 people. The term core-based statistical area (CBSA) refers to both these micro and metro areas. CBSAs are made of counties, including the county containing the urban core and any adjacent counties with “a high degree of social and economic integration” (measured by commute to work).1
Use the core_based_statistical_areas() function to get boundary data for all metropolitan and micropolitan areas in the United States. To find a specific place, filter the NAME column; in this case we are looking for Memphis.
metro <- core_based_statistical_areas() %>%
filter(str_detect(NAME, "Memphis"))
#' same thing in base R notation:
#' cbsa <- core_based_statistical_areas()
#' metro <- cbsa[grep("Memphis", cbsa$NAME), ]
plot(metro$geometry)
Counties
The Memphis metro area covers three states: Tennessee, Mississippi, and Arkansas. If we want to see which counties are in the metro area, we’ll need to get data for each state. We can do with the counties() function.
To get data for multiple states at once, we can combine the map_dfr() and counties functions. This tells R to get data from "AR", "TN", and "MS", then joins those tables into one: stateCounties.
The result, stateCounties (not shown) is a map of all counties in Arkansas, Tennessee, and Mississippi. We can then use base R notation to filter to only counties within the metro area (sf package required).
metroCounties <- stateCounties[metro, op = st_within]
#' Note: op = st_within only returns counties within the metro area.
#' Use stateCounties[metro, ] to return all counties that are within or border the metro area.
To create a map with multiple customized layers, I’ll use ggplot. Below, the Memphis metro area is outlined in red and metroCounties are outlined in grey.
ggplot() +
geom_sf(data = metroCounties, fill = "white", color = "grey") +
geom_sf(data = metro, fill = NA, color = "red") +
theme_void()
Figure 1: Counties in Memphis Metro Area
Places
Similarly, we can all the cities, towns, and villages within the metro area. A place is either a legally incorporated area or a Census Designated Place (CDP), which is not legally incorporated and is used for statistics purposes. There’s no population size requirements for CDPs.
Using R we can easily copy and paste the above process and make small adjustments to create a new plot for all the places within the Memphis metro area.
statePlaces <- map_dfr(c("AR", "TN", "MS"), ~{
places(.x)
})
metroPlaces <- statePlaces[metro, op = st_within]
ggplot() +
geom_sf(data = metroPlaces, fill = "white", color = "grey") +
geom_sf(data = metro, fill = NA, color = "red") +
theme_void()
Figure 2: Places in Memphis metro area
Interactive Mapview Map
Making static maps is easy enough, but R packages also makes it simple to create interactive maps. For this first map I’ll be using the mapview package. mapview is great for creating simple maps with basic useful features.
It’s possible to create a map with just mapview(metro), and adding another layer is as simple as + mapview(metroCounties). I went ahead and adjusted some of the aesthetics for the below map.
Below is a basic interactive map of the Memphis Metro Area, including counties and places within the metro area.
mapview(metro,
layer.name = "Memphis Metro Area",
label = metro$NAMELSAD,
col.regions = "#EC9061", #' fill color
alpha.regions = 0.1, #' fill opacity
lwd = 3, #' line width
color = "#EC9061" #' line color
) +
mapview(metroCounties,
layer.name = "Counties in Metro",
label = metroCounties$NAMELSAD,
col.regions = "#F8F5E6",
alpha.regions = 0.3
) +
mapview(metroPlaces,
layer.name = "Places in Metro",
label = metroPlaces$NAME,
col.regions = "#46ABBF",
alpha.regions = 0.4
)
The layers button in the upper left, default in mapview, allows you to change the base layer or turn layers on and off. Hovering over a region displays the Census label name, set in the above code by label = dataSet$variable. mapview also enables popups by default; clicking on a feature displays all available attribute data.
Urban Areas
CBSA boundaries are set by county lines rather than where people actually live and work, meaning their boundaries can include too much rural land for what we’re trying to study. Urban areas include “urbanized areas,” which are densely developed areas with at least 50,000 people, and “urban clusters,” which have a between 2,500 and 50,000 people.2
urb <- urban_areas() %>%
filter(str_detect(NAME10, "Memphis"))
#' alt example in base R notation:
#' uas <- urban_areas()
#' urb <- uas[grep("Memphis", uas$NAME10), ]
ggplot() +
geom_sf(data = metroCounties, fill = "white", color = "grey") +
geom_sf(data = urb, color = "red") +
theme_void()
Figure 3: Memphis urban area
After overlaying the maps, we can see the Memphis urban area is significantly smaller than the overall metro boundaries. Areas outside the boundary are considered rural. We can also see gaps within the urban area where no people would live, like the airport and Shelby Farms.
ZCTAs
Zip Code Tabulation Areas (ZCTAs) are Census representations of USPS ZIP Codes.3 They are created using Census blocks and are not bound by place, county, tract, or block group.
Use zcta() to download all ZCTAs in the US, filterable by state. Because ZIP Codes starting with the same two digits will usually be clustered together, you can also filter by starts_with(). The ZCTA file is massive, so it’s a good idea to filter one way or another.
While I could use map_dfr() again, the Memphis urban area is actually used in the CRAN example for how to use the zctas() function. According to the example, ZCTAs in the Memphis area start with “37,” “38,” and “72.” After working through the example, I found I only needed “38” and “72” (“37” covered Middle Tennessee). Then I further filtered to only keep ZCTAs intersecting the Memphis urban area.
zcta <- zctas(cb = TRUE, starts_with = c("38", "72"))
urbZCTA <- zcta[urb, ]
ggplot() +
geom_sf(data = metroCounties, fill = "white", color = "grey") +
geom_sf(data = urbZCTA, fill = NA, color = "red") +
theme_void()
Figure 4: ZCTAs in Memphis urban area
Notice how ZCTAs do not always align with county boundaries.
PUMAs
Public Use Microdata Areas (PUMAs) are Census boundaries set every 10 years which contain at least 100,000 people. They are notably used with Public Use Microdata Sample (PUMS) data, which are anonymized individual-level Census records.4 This data is useful to researchers who want to create custom queries of data rather than using pre-tabulated estimates provided by the Census Bureau.
statePUMAS <- map_dfr(c("AR", "TN", "MS"), ~{
pumas(.x)
})
urbPUMAS <- statePUMAS[urb, ]
ggplot() +
geom_sf(data = metroCounties, fill = "white", color = "grey") +
geom_sf(data = urbPUMAS, fill = NA, color = "red") +
theme_void()
Figure 5: PUMAs in Memphis urban area
PUMAs boundaries are built on census tracts. Though they may be guided by county lines, they are not limited to them, and one PUMA can cover multiple counties.
In the above map, I was concerned that there was no PUMA touching the southeast corner of Shelby County (center). I ran the code again, instead filtering by metroCounties, so I could see additional PUMAs.
metroPUMAS <- statePUMAS[metroCounties, ]
plot(metroPUMAS$geometry)
Figure 6: PUMAs in Memphis metro area
The PUMA touching the SE of Shelby County extends the remainder of north Mississippi, far outside the area I’m trying to study. So, if using PUMAs to study the Memphis metro, I would use urbPUMAS.
Interactive Leaflet Map
Originally I tried making the plots in this section into another interactive map with mapview, like in ??. However, the boundaries in this section do not fit neatly into one another. When stacked, they look messy. Users could turn off layers in the settings, but I think this is expecting a lot out of users; if the map is initially busy/confusing, they won’t want to use it at all.
At this point, I decided to switch to the leaflet package, which is easier for maps requiring more customizing.
Leaflet’s syntax takes a little bit more to get started. Instead of mapview() + mapview(), the minimum code for two layers is leaflet() %>% addTiles() %>% addPolygons() %>% addPolygons().
To show and hide layers, first assign a group to each feature. Enable addLayersControl() and list which groups you want users to be able to toggle on/off. To hide specific groups initially, use hideGroup(). Read more in the Leaflet for R guide.
The following is an interactive map of the Memphis urban area, with surrounding ZCTAs and PUMAs.
Show code
urbMap <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB (default)") %>%
addPolygons(data = urb,
fillColor = "#c93335",
fillOpacity = 0.2,
weight = 1,
color = "#c93335",
opacity = 1,
label = ~NAME10,
group = "Urban Area") %>%
addPolygons(data = urbPUMAS,
fillColor = "#9379c2",
fillOpacity = 0.2,
weight = 1,
color = "#9379c2",
opacity = 1,
label = ~NAMELSAD10,
group = "PUMAs",
highlightOptions = highlightOptions(
color = "white",
weight = 3,
bringToFront = TRUE)) %>%
addPolygons(data = urbZCTA,
fillColor = "#66b888",
fillOpacity = 0.2,
weight = 1,
color = "#66b888",
opacity = 1,
label = ~ZCTA5CE10,
group = "ZCTAs",
highlightOptions = highlightOptions(
color = "white",
weight = 3,
bringToFront = TRUE)) %>%
addLayersControl(overlayGroups = c("Urban Area", "PUMAs", "ZCTAs"),
options = layersControlOptions(collapsed = FALSE)) %>%
addScaleBar(position = "bottomleft")
urbMap %>%
hideGroup(c("ZCTAs", "PUMAs"))
Figure 7: Interactive leaflet map of Memphis urban area
Shelby County
Next let’s focus on Shelby County, which contains Memphis. We can filter our metroCounties data from earlier.
shelby <- metroCounties %>%
filter(str_detect(NAME, "Shelby"))
This gives us an outline of Shelby County (not shown here).
Tracts and Block Groups
We can get Census tract data using the tracts() function. At this geographic level we can also specify a county name, such as “Shelby,” so we don’t have to filter the data.
tracts <- tracts("TN", "Shelby")
blkgrp <- block_groups("TN", "Shelby")
ggplot() +
geom_sf(data = tracts, fill = "white", color = "red") +
geom_sf(data = shelby, fill = NA, color = "black") +
theme_void()
ggplot() +
geom_sf(data = blkgrp, fill = "white", color = "grey") +
geom_sf(data = tracts, fill = NA, color = "red") +
geom_sf(data = shelby, fill = NA, color = "black") +
theme_void()
Figure 8: Census tracts and block groups in Shelby County
School Districts
School districts have been a contentious topic over the past decade, after Memphis City Schools dissolved its charter to merge into the county school district. Following this, suburban towns like Collierville and Germantown separated from Shelby County Schools and formed their own municipal school districts.
The following code gets all school districts in Shelby County
stateSchools <- school_districts("TN")
schools <- stateSchools[shelby, op = st_within]
Final Map
Below is an interactive map combining all of the maps created in this post.
Show code
final <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "CartoDB (default)") %>%
addTiles(group = "OSM") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "ESRI World Imagery") %>%
addPolygons(
data = metro,
group = "Metro Area",
label = metro$NAMELSAD,
fillColor = "#39383d",
fillOpacity = 0.2,
weight = 4,
opacity = 1,
color = "#39383d",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = metroCounties,
group = "Counties",
label = metroCounties$NAMELSAD,
fillColor = "#c93335",
fillOpacity = 0.2,
weight = 3,
opacity = 1,
color = "#c93335",
popup = TRUE,
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = metroPlaces,
group = "Places",
label = metroPlaces$NAME,
fillColor = "#9379c2",
fillOpacity = 0.2,
weight = 2,
opacity = 1,
color = "#9379c2",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = urb,
group = "Urban Area",
label = urb$NAME10,
fillColor = "#71b6f7",
fillOpacity = 0.2,
weight = 3,
opacity = 1,
color = "#71b6f7",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = urbZCTA,
group = "ZCTAs",
label = urbZCTA$ZCTA5CE10,
fillColor = "#66b888",
fillOpacity = 0.2,
weight = 3,
opacity = 1,
color = "#66b888",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = urbPUMAS,
group = "PUMAs",
label = urbPUMAS$NAMELSAD10,
fillColor = "#a4d740",
fillOpacity = 0.2,
weight = 2,
opacity = 1,
color = "#a4d740",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = tracts,
group = "Census Tracts",
label = tracts$NAMELSAD,
fillColor = "#ee6668",
fillOpacity = 0.2,
weight = 2.5,
opacity = 1,
color = "#ee6668",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = blkgrp,
group = "Block Groups",
label = blkgrp$GEOID,
fillColor = "#f28b73",
fillOpacity = 0.2,
weight = 1,
opacity = 1,
color = "#f28b73",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addPolygons(
data = schools,
group = "School Districts",
label = schools$NAME,
fillColor = "#dbb4f8",
fillOpacity = 0.2,
weight = 2,
opacity = 1,
color = "#dbb4f8",
highlightOptions = highlightOptions(fillOpacity = 0.5)
) %>%
addLayersControl(
baseGroups = c("CartoDB (default)", "OSM", "ESRI World Imagery"),
overlayGroups = c(
"Metro Area",
"Counties",
"Census Tracts",
"Block Groups",
"ZCTAs",
"PUMAs",
"Urban Area",
"Places",
"School Districts"
),
options = layersControlOptions(collapsed = FALSE)
) %>%
addScaleBar(position = "bottomleft")
final %>% hideGroup(
c(
"Metro Area",
"Places",
"Urban Area",
"ZCTAs",
"PUMAs",
"Census Tracts",
"Block Groups",
"School Districts"
)
)
Figure 9: Interactive leaflet map of Memphis boundaries
Acknowledgments
To create this post I mostly referenced the book Analyzing US Census Data by Kyle Walker, particularly Chapter 5 and the Chapter 7 section on spatial overlay.
The leaflet website and package documentation were useful when I got stuck.
Color palette for maps were created using this Animal Crossing: New Horizons palette by sokea-cc and by referencing a screenshot of the NookPhone from AC:NH on Nookipedia.
US Census Bureau, “ZIP Code Tabulation Areas (ZCTAs),” October 8, 2021, https://www.census.gov/programs-surveys/geography/guidance/geo-areas/zctas.html.↩︎
Tigris, “Urban_areas: Download an Urban Areas Shapefile into R in Tigris: Load Census TIGER/Line Shapefiles,” September 23, 2021, https://rdrr.io/cran/tigris/man/urban_areas.html.↩︎
Until the API is complete, PUMAs data is available for download from https://www.ipums.org/.↩︎